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O . 

^S) , We study the dynamics of bright and dark matter-wave solitons in the presence of a spatially 

varying nonlinearity. When the spatial variation does not involve zero crossings, a transformation 
is used to bring the problem to a standard nonlinear Schrodinger form, but with two additional 
) ^ terms: an effective potential one and a non-potential term. We illustrate how to apply perturbation 

, theory of dark and bright solitons to the transformed equations. We develop the general case, but 

' primarily focus on the non-standard special case whereby the potential term vanishes, for an inverse 

square spatial dependence of the nonlinearity. In both cases of repulsive and attractive interactions, 
appropriate versions of the soliton perturbation theory are shown to accurately describe the soliton 
dynamics. 

Ph! PACS numbers: 03.75.Lm, 05.45.Yv, 02.30.Jr 

O^. I. INTRODUCTION 

The experimental creation of atomic Bose-Einstein condensates (BECs) has been one of the most 
^ \ fundamental developments in quantum and atomic physics over the past two decades. The impressive 

. progress in this field due to intense experimental and theoretical studies has been already summarized 

I in various books P, H, [l] and reviews [1]. This progress has been, to a considerable extent, fueled 

by the fact that, in a mean-field picture, BECs can be described by a macroscopic wavefunction 
obeying the Gross-Pitaevskii (GP) equation, which is an equation of the nonlinear Schrodinger 
(NLS) type. In such a mean-field description, the effective nonlinearity (which is introduced by 
■ interatomic interactions) allows for studies of macroscopic nonlinear matter waves; in this respect 

\ it is important to note that bright matter- wave solitons in attractive BECs 0, 0] , as well as dark 

• • ■ d, H, [13, [O, [H, [H, 13 aiid gap il3 matter wave solitons in repulsive BECs, have been observed in 

_ ^ a series of experiments (see also the recent review Q ) • 

. One of the remarkable possibilities arising in the physics of BECs is that the interatomic inter- 

' actions (and, hence, the effective nonlinearity) can be manipulated by means of different types of 

5^ , temporally- or spatially- varying external potentials. More specifically, the s-wave scattering length 

(which is proportional to the nonlinearity coefficient in the GP equation) can be experimentally ad- 
justed using either magnetic [l^, [l3| or optical Feshbach resonances [l^ in a very broad range. The 
availability of these tools has led to a number of consecutive theoretical and experimental studies. 
For instance, the formation of bright matter-wave solitons and soliton trains of ^Li 0, Q and ®^Rb 
Q atoms used a tuning of the interatomic interactions from repulsive to attractive. Also, this type 
of manipulations was instrumental in achieving the formation of molecular condensates |19| , and the 
probing of the BEC-BCS crossover A parallel track of theoretical studies has explored the use of 
a time-dependent modulation of the nonlinearity coefficient to stabilize attractive higher-dimensional 
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BECs against collapse 2l|, or to create robust matter- wave breathers in lower-dimensional BECs 



[22( 1 . More recently, the use of spatial variations of the nonlinearity to create so-called "collisionally 
inhomogeneous" environments has been pro posed. In that regard, major developments included adi- 
abatic compression of matter- waves (23l. |24| . Bloch oscillations of matter- wave solitons [23*1, atomic 
soliton emission and atom lasers [25| . enhancement of transmittivity of matter- waves through bar- 
riers [2^ [l^] , dynamical trapping of matter- wave solitons [l^ , stable condensates exhibiting both 
attractive and repulsive interatomic interactions [28l| as well as the delocalization transition of mat- 
ter waves [29i] . Among the types of spatial variations of the nonlinearity that have been proposed, 
one can trace linear one s [23. [26ll . a s well as parabolic [30| , random [Sll , periodic [1^, [H, [H, H^l , 
and localized (step-like) [25l l35l. l3a| ones. On the mathematical side, a number of detailed studies 
[37I m, [3^ have appeared, addressing aspects such as the effect of a "nonlinear lattice potential" 
(i.e., a spatially periodic nonlinearity) on the stability of matter-wave solitons, and the interplay 
between drift and diffraction/blow-up instabilities. More recently, the interplay of nonlinear and 
linear potentials has been examined in both continuum [i^l and discrete [4l| settings (see also the 
recent work ^i^l and references therein). 

Our aim in this work is to study the dynamics of matter- wave solitons in the presence of a spatially- 
dependent nonlinearity. We consider both dark solitons in repulsive BECs, as well as bright solitons 
in attractive BECs. In the case where the sign of the nonlinearity coefficient (hereafter referred 
to as g{x)) does not change, we first show that a change of variables can convert the spatially 
variable nonlinearity problem into a "regular" one where the nonlinearity has a constant prefactor. 
This transformation results in the emergence of two additional perturbation terms: one of them 
can be considered as an effective potential term (i.e., a spatially-dependent function multiplying 
the macroscopic wavefunction u), while the other one can not (it consists of a spatially-dependent 
function multiplying the derivative of the wavefunction dxU, for an elongated BEC along the x- 
direction). We use this transformation as a starting point in order to develop perturbation theory 
for the soliton dynamics in the presence of g(x) for the case of arbitrary g{x). However our focus 
is on the case where g(x) is such that the potential term completely vanishes. The reason for this 
selection is that it appears to be the less physically intuitive case (due to the derivative nature of 
the corresponding perturbation). Moreover, the effect of a perturbation induced by a "standard" 
potential term has been studied fairly extensively in the BEC context (see, e.g., [ESSBDj also in 
the particular case of coUisionally inhomogeneous BECs (see, e.g., [H, [H, [26[ ) . 

Our investigation is structured as follows. In the next section, we give the general setting and 
analyze the relevant transformation. In section III, we focus on dark solitons, first providing the 
general theory, and then applying it to the particular case of interest. In section IV, we follow a 
similar path for the case of bright solitons. Finally, in section V, we summarize our findings and 
present our conclusions, as well as some interesting directions for future study. 



II. PERTURBED GROSS-PITAEVSKII EQUATION AND THE DERIVATIVE-ONLY 

CASE 

In this work we will restrict ourselves to an effective one-dimensional (ID) description accounting 
only for the longitudinal dynamics of the condensate. In the transverse directions the atoms should 
be tightly confined which can be realized by an isotropic harmonic potential with a trap frequency 
L0± (associated with the harmonic oscillator length a±). For a± small enough, one can regard the 
transversal dynamics as frozen (see chap. 1 in Q and also [H, SI]), i.e., only the corresponding 
ground state is occupied. The longitudinal motion takes place in the x-direction and should not be 
confined. Then, the corresponding mean-field equation is given by 

^2 

ihdt-^ ^ -—dl'^ + g{x)\'^f'^ (1) 
2m 

where '^{x,t) is the macroscopic wave function, m the atomic mass and g{x) = g^^^^ {x)/2TTa'^ 
the effective ID interaction coefficient. The parameter g^^^\x) — Airh^a/m characterizes the two- 
particle interaction in 3D with the s-wave scattering length a. The latter is positive (negative) for 
repulsive (attractive) condensates consisting of, e.g., ^^Na (''Li) atoms. The value of the scattering 
length can be tuned, as mentioned above, e.g., by use of magnetic Feshbach resonances [i^. In the 
vicinity of a magnetic Feshbach resonance the value of the scattering length depends on the value of 
an applied magnetic field. Thus, one can achieve a spatially dependent scattering length by applying 
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an inhomogeneous magnetic field yielding a coUisionally inhomogeneous BEC. We now use suitable 
straightforward rescalings (see, e.g., [2^) and dimcnsionlcss units to express Eq. ^ in the following 
form: 

zdt^^-^dl^ + s\gix)mH, (2) 

where the coefficient s = sgn((7) = ±1 for attractive and repulsive condensates, respectively. Apply- 
ing the transformation = allows us to rewrite eq. ([2]) in the following way: 

idtu = + s\ufu + Veff{x)u - ^/gdx — dxU, (3) 

with the effective potential term Veff{x) = ~\^d^^- Equation ([3]) can be written as the usual 
NLS equation (with a defocusing or focusing nonlinearity for s = ±1, respectively) with an external 
spatially-dependent perturbation P[u{x,t)]x\, namely: 

idtU + ^d^u - s|upu = P[u{x,t);x]. (4) 

The perturbation can be expressed as P[u{x, t)]x] — Pi^[u{x, t); x] +Ppfp[u(x, t);x], i.e., it consists of 
a linear effective potential contribution PL[u{x,t); x] = Veff{x)u{x,t), as well as of a non-potential 
perturbation of the form Pi\fp[u{x,t);x] = —y/gdx^dxU. In this work we are mainly interested 

in the effects of the less standard, non-potential type of perturbation. Therefore, we assume the 
coUisionally inhomogeneous interaction to be of the form 

^ (D+cxr ^ 

with arbitrary constants C and D. For such a particular selection oi g{x), Veff{x) vanishes leading 
to the perturbation 

C 

P[u{x,t);x] = -———dxU, (6) 
U + Cx 

consisting only of the non-potential contribution in the right hand side. Thus we can investigate the 
pure effects of the latter, non-standard contribution (the effects of a standard linear potential have 
been studied fairly extensively; see e.g. [H, 0, H, 01 ) ■ We choose C — 1 and D = —200; thereby the 
singularity in the perturbation occurs at ccg*"® = 200 (which will be outside the region of interest in 
our domain). For this choice, the spatial dependence of the coefficient g{x) is shown in Fig. [TJ 
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FIG. 1: Spatial dependence of the interaction parameter g{x) for C — 1 and D = —200 for a repulsive 
condensate. 



A Full perturbative approach 



4 



III. DARK MATTER- WAVE SOLITONS 
A. Full perturbative approach 

Let us first consider the case of dark matter- wave solitons for s = 1. In order to treat effects of the 
perturbation ([H]) on a dark soliton analytically we employ the adiabatic perturbation theory assuming 
that the functional form of the soliton remains unchanged by the perturbation (an assumption that 
in our setting will be justified a posteriori). We first use the transformation u —>■ uexp{—it) to put 
Eq. (HI) in the form idtu + ^d^u — — l)u — P[u{x,t);x] and use as an ansatz for the soliton, 
the following expression, 

u = BtSinh{B{x - Xq)) + lA, (7) 

which is the exact dark soliton solution of the above mentioned unperturbed NLS equation. Accord- 
ing to the above discussion, the soliton depth A, velocity B (with A"^ + B'^ = 1), and center xq are 
assumed to be unknown functions of time. The Lagrangian density of an unperturbed dark soliton 
is given by [isj : 

du) = '-iu*dtu - udtu*){l l\d.uf - ^i\u\' - l)^ (8) 

while the averaged Lagrangian, L — J dx£{u), can be calculated by substituting ([7]) in Eq. ([8]) 
yielding: 

LiA, xo) = 2dtXo{-AB + tan-i(^)) - ^B^ (9) 

In (47j it was shown that, within the framework of the adiabatic perturbation theory for small 
perturbations, the parameters of the soliton Uj = {xo,A} obey the the Euler-Lagrange equations 

da L - -^da,' L = 2Re{ f dxP*{u)da .u} (10) 

with a'j — dtaj. This leads to a system of ordinary differential equations (ODE) for A and xq: 

dtA = ^B^ J dxy^d,^sech'^{B{x - xo)) 

+ -B^ / dxy^dl—taiih{B{x - xo))sech^{B{x - Xo)), (11) 
4 J v5 

c'f.To = A—^-A j dx^/gdx—psech^{B{x~xo))(ta,nh{B{x~xo)) 
+B(x — a;o)sech^(i3(a; — a^o))^ 
-i y dx^gdl^Jy{^iKv.\,\B{x-xo))-\) 

■^{x - Xo) tanh(B(a:: - Xo))^&(A? [B{x - Xq))^ . (12) 

For an interaction obeying Eq. ([5]), the terms arising from the linear potential vanish, leading to 
the system: 

dtA - f dx—^sech^{B{x-xo)), (13) 

2 J D + Cx 

1 f C / 

dtXo = A— -A / dx— — — sech^(B(a; — a;o)) ( tanh(i3(a; — Xo)) 

2 J D + Cx V 

+B{x-xo)sech^{B{x-xo))y (14) 
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B. Approximations 



We can simplify the general framework of Eqs. (|llll2p by performing a Taylor expansion of the 
interaction around x = xq^ leading in first order to 

1 



dtA « (1-^2)7^9,—== (15) 
I ^ I 

dtxo w A+ -—^.yg{xo)dl—-^ . (16) 

4 iJ^ g[x) x=xo 

The Taylor expansion around xo can be justified in most settings due to the exponential localization 
of the soliton around its center. Dropping higher-order terms essentially implies that the interaction 
does not change on the scale of the width of the soliton (if the latter assumption is invalid, then 
we can not resort to this approximation). In the special case of the interaction ([5]) the following 
evolution equations are obtained: 

dtXo « A (18) 

By combination of Eqs. (fT7|) and (HH]), we obtain a single second-order ordinary differential equation 
(ODE) 

for the center of the soliton. If the velocity of the soliton is small, dtXQ <C 1, one can neglect the 
second term in the right hand side of Eq. leading to: 



We will discuss the validity of this approximation in our numerical results below. Equation (I20|l is 
the equation of motion (EOM) 

d^xo = ~d^„V^"{xo), (21) 
of a particle in the presence of the effective potential 

V^"{xo) = -^ln(|C.To + Z?|). (22) 

Therefore, we will denote Eq. as EOM and Eq. pp)) as EOMa in the next section. As an 

interesting aside, we should note that even in the presence of the kinetic term (i.e., if (dtXo)^ is not 
neglected), one can rewrite Eq. (fT9|) as a Hamiltonian system (see Ref. [11]) using a generalized 
momentum P = g{xo)dtXo. With this momentum one finds for a system 

d^xo = fixo)il-idtXof), (23) 

the equations of motion 

dtXo = (24) 

d,P = 3MfM + P{j^-^), (25) 



which correspond to the Hamiltonian 



Hixo,P)=^l^^+F{xo), (26) 
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with 

/Xo 
dx'ofix'o)), (27) 

Fixo) = -Igixo). (28) 

In the particular case of f{xo) — f d+Cxo ' ^^"^ obtains for the momentum 

P = AiD + Cxo)^dtXo, (29) 

and for the Hamiltonian 

Hixo,P) = QiD + Cxo)-^ - ^AiD + Cxo)-^. (30) 

C. Numerical Results 

In this section we present and compare the numerical results obtained by solving the full partial 
differential equation (PDE) of the GP type as weU as the ODEs p3ll4p . the EOM ^ and the 
simplified EOMa (|20p . We have confirmed that throughout our simulations the soliton is localized in 
a region with a well defined perturbation avoiding the singularity. Since the soliton is exponentially 
localized, the spatial integrations in Eqs. (|I3|I4p can be restricted to a region around the center of 
the soliton with a finite perturbation and a well defined integrand. The time evolution is performed 
by the Adams-Bashforth-Moulton predictor-corrector method. 

Figure[2]shows the time evolution of the density profile of a dark soliton with xo(0) = for different 
initial velocities obtained by solving Eq. (jl]). Black represents the highest density, while white 
corresponds to the lowest density. The dotted lines are the results for xq (the center of the soliton) 
as obtained by solving the ODEs (|13|14p . The results agree very well, showing that the adiabatic 
perturbation theory describes the motion of the center of the soliton accurately. For A^nit — 
(lowest curve) the soliton gets accelerated to the negative half-plane and moves immediately into 
this direction. For Amit = 0.25 and Ainu — 0.5 (middle and top curve, respectively), the soliton 
also gets accelerated into the direction of the negative half-plane but starts moving to the positive 
one due to its initial velocity until it reaches a turning point of zero velocity and changes direction. 
By investigation of Fig. [T] one observes that in the considered region the value of the interaction 
parameter decreases for decreasing x. Thus, the soliton gets accelerated into the direction with 
a smaller interaction parameter. Due to the fact that the interaction is repulsive, the interaction 
energy decreases with decreasing interaction parameter. So the soliton tends to move into the region 
with less interaction energy. This happens despite the fact that the interaction parameter does 
not enter explicitly in the perturbation as a potential but rather through the product of its first 
derivative and the first derivative of the scaled wavefunction. 




20 40 60 80 100 
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FIG. 2: Time evolution of the density of dark solitons with a;o(0) — and Ainu — 0, 0.25, 0.5 (from bottom 
to top). The dotted hne is the CM parameter (xo) obtained by solving the ODEs (I13p - (|14| l. 
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Figure [3] shows the time evohition for xq{0) = 0, Ainu = 0(a) and Ainn = 0.5(b) of the differences 
between the results for the center of the sohton obtained by solving the PDE, and the ODE, EOM 
and EOMa, respectively. We calculated the center of mass of the PDE solution by performing the 
integration xq = J x{h — \u\'^)dx/ J{b — \u\'^)dx, with b = \u{xb)\'^ being the background density 
evaluated far away from the center of the soliton. The differences for the ODEs and the EOM 
are almost equal and small for both initial velocities. So the adiabatic perturbation theory works 
fine for describing the center of the soliton. The results obtained by solving EOMq coincide for a 
small time period with the result of the PDE. For longer times, they deviate from these results. 
For a larger initial velocity the deviation is even larger. The reason for this is that we neglected 
the impact of the velocity in EOM^ and, thus, the approximation gets worse for larger velocities. 
However, the qualitative behavior is described correctly even within this approximation. Hence, one 
can understand the behavior of the soliton as a particle moving in the effective potential (|22p and 
thus can explain the acceleration observed in Fig. [2l 




100 



t t 

FIG. 3: Difference of the soliton center Aa;o calculated by solving the PDE with the result of the ODEs, 
EOM and EOM„ for Ai^u = 0.0 (a) and Ai^u = 0.5 (b). 



IV. BRIGHT MATTER- WAVE SOLITONS 



A. Full perturbative approach 

In the case of attractive interactions (s — —1), Eq. (|3|) reads after substitution of t = t/2 

idrU + dlu + 2\u\^u = 2eP{u). (31) 

In the absence of perturbations, it is well known that Eq. (pij) possesses a bright soliton solution of 
the form 

u{z,t) = 2i77exp(-2i^a; - i$)sech(z), (32) 

where z — 2ri{x—Q, while rj represents the amplitude, the phase and C the center of the soliton, and 
^ is related to the velocity of the soliton. For a small perturbation, we can now employ the adiabatic 
perturbation theory for bright solitons [49l| to treat the perturbation effects analytically. Then, 
the soliton parameters become slowly-varying functions of time, however the shape of the soliton 
remains unchanged (once again this is a principal assumption that will be justified a posteriori). 
With the general perturbation arising due to a spatially-dependent scattering length, one arrives at 
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the following system of ordinary differential equations for the parameters of the soliton: 

drV = 877^^ J dx^d^^sech^{2Tj{x - 0), (33) 

J V 5 

-27f [ dx^dl^ta.nh(2r](x - C)sech'^(2r](x - ()), (34) 
J V9 

drC = -4^ + 877$ J dx^d^^ix - Osech^i27^{x - 0), (35) 



4(C^ - V^) + 8?7^ / dxy^d^^sech^(2'n(x - 0) tanh(2?7(a; - C)) 

J Va 

(1 - 2rjxta,nh{2ri{x ~ C))) 

-2r/ f dx^dl—sech'^{2f]{x - C))(l - 2?72;tanh(2r/(a; - C))). (36) 
J \f9 

For an interaction parameter of the form (O one obtains: 

drn = 877^^ 1 d:,_^sech2(2r;(x~C)) (37) 

dri = 87]^ / dx ^ tanh^(277(a:^C))sech^(27?(a;-C)) (38) 

J U + Cx 

drC = -4? + 877^ j dx^f^ix - C)sech2(27y(x - 0) (39) 

= 4(^2 _ ^ /■ ^L_sech^ (277(2; - C)) tanh(2r/(a; - C)) 

J Z) + ox 

(1 - 2?7a; tanh(277(a; -()))■ (40) 
B. Approximations 

From the above equations it is clear that Eq. (f40|) describes the time evolution of the phase of 
the soliton which, however, does not emerge in the equations determining the other parameters. 
Therefore, we will restrict our considerations to Eqs. ([37| - p9p . Since the soliton is exponentially 
localized around x — we can perform a Taylor expansion around C and thus simplify Eqs. (|33|) -(|35 p 
as follows: 

1 



drV = 87J^yMOd,^= ^ (41) 



1 



3' 



dr^ = -V^yMOd.^= ^ (42) 



x=C 



drC = -4e (43) 

The physical interpretation of this approximation is that the interaction parameter does not vary 
over the width of the soliton. Focusing more specifically on a collisional inhomogeneity of the form 
of Eq. lU where these contributions vanish exactly, leads to 

3tC = -4{. (46) 
We can solve the simplified Eq. (|44|l directly by using Eq. (|46|) : 



(47) 
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Combination of Eqs. (|45ll47p and back transformation to the real time t leads to the equation of 
motion for the soliton center: 

with an associated effective potential: 



(49) 



C. Numerical Results 



Figure m shows the time evolution of the density of the bright soliton (note that, in contrary to 
before, black represent zero density, while white represents high density). The results are obtained 
by integration of Eq. ([4]) with an initial state given by Eq. (j32p with parameters initialized as 
T] — 0.5, <i> = C = and (,init = 0, 0.25, 0.5 (from top to bottom). The dotted line shows the 
corresponding results for the center of the soliton ( obtained by solving Eqs. l|37M0p . The results 
of the perturbation theory once again agree very well with the results of the PDE. In the case of 
zero initial velocity (top curve), the soliton gets accelerated to the positive half- plane and starts 
moving into this direction immediately. For positive S,init the initial velocity is negative, according 
to Eq. ((38|) . leading to a motion towards the negative half-plane. However, the soliton still moves 
toward the direction of the positive half-plane, due to its initial speed, yet eventually it acquires a 
zero velocity and a change of the direction of motion occurs. The direction of the acceleration is the 
direction of increasing interaction parameter as can be seen by comparing the results with Fig. [H as 
this minimizes the energy of the system (even though the interaction does not act, strictly speaking, 
as a potential). 
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FIG. 4: Time evolution of the density of a bright soliton with initial parameter rj — 0.5, $ — = and 
£,init = 0, 0.25, 0.5 (from top to bottom). The dotted line shows the result for ( (the center of the soliton) 
from the adiabatic perturbation theory. 

Figure [5] shows the difference of the center of the soliton calculated by the PDE with the results of 
the ODEs dSZIl-lISHl) and the EOM (gS]). The center of the soliton of the PDE solution is determined 
by the quotient C = / x\u\'^dx/ J \u\^dx. Figure[5^ shows the differences for 77 = 0.5, $ = C = and 
^init — 0. The difference of the EOM result from the PDE result is slightly larger than the difference 
of the ODE result. Both differences increase with time but they are still very small for the time 
period considered. Fig. shows the differences for 77 — 0.5, <f> = C = and ^i„it — 0.5. In this 
case, the absolute differences are an order of magnitude larger than in the previous case. However, 
compared to the position and width of the soliton one can still regard them as small. The results of 
the ODEs and the EOM are almost equal. A conclusion of this investigation is that the dynamics of 
the soliton is described fairly accurately by the model of a particle subject to the effective potential 
dMl). 
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FIG. 5: Difference of the solitori center calculated by solving the ODE and EOM to the result of the PDE 
with initial parameter r; = 0.5, $ = ^ = and ^init = (a) or ^init = 0.5 (b). 



V. CONCLUSIONS AND FUTURE CHALLENGES 



In this work, wc considered the effect of (slowly-varying) spatially dependent nonlinearities of a 
definite sign on both dark and bright matter- wave solitons of repulsive and attractive Bose-Einstein 
condensates, respectively. Wc have shown that a relevant transformation can be employed to con- 
vert the spatially dependent nonlinear problem into one of spatially uniform nonlinearity, at the 
expense of introducing two perturbative terms. One of the latter is in the form of a linear potential 
(which have been considered extensively previously), while the other constitutes a non-potential 
type of perturbation, being proportional to the spatial derivative of the field. To especially highlight 
the non-potential nature of the second term, we considered collisional inhomogeneities of inverse 
square spatial dependence, whereby the linear potential perturbation identically vanishes, and the 
purely non-potential one has to be considered. Even in these settings (but also more generally), we 
found that soliton perturbation theory provides a powerful tool towards describing such collisional 
inhomogeneities . 

It would be interesting to extend the present considerations in a number of directions. Firstly, it 
would be relevant to appreciate the effect of g{x) on higher-dimensional structures, such as vortices, 
and on their stability. On the other hand, it would be especially interesting even in one spatial 
dimension to determine whether techniques like the ones used here (or variants thereof) can be 
applied to cases where the sign of the nonlinearity changes. Finally, it woidd be relevant to observe 
systematically how techniques such as soliton perturbation theory may fail, as the size of spatial 
extent of the collisional inhomogeneity becomes comparable to that of the solitary wave and to 
understand the ensuing phenomenology in such cases. Studies along some of these directions are 
currently in progress and will be reported in future publications. 
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